#include "hls_img.h"

Hls_img::Hls_img()
{
    //ctor
}

Hls_img::~Hls_img()
{
    //dtor
}

//--------------------------------------------------------------------------
//Statistical Geometrical Features
//Features Function

//--------------------------------------------------------------------------
int Hls_img::isSafe(Mat M, int row, int col, Mat visited, bool x, Mat matric, int count)
{
    Scalar intensity1 = M.at<uchar>(row, col);
    Scalar intensity2 = visited.at<uchar>(row, col);
    if((row >= 0) && (row < M.rows) && (col >= 0) && (col < M.cols)
       && (intensity1.val[0] == x && intensity2.val[0] == 0))
    {
        matric.at<uchar>(row, col) = count;
        return 1;
    }
    return 0;
}

void Hls_img::DFS(Mat M, int row, int col, Mat visited, bool x, Mat matric, int count)
{
    // These arrays are used to get row and column numbers of 8 neighbors
    // of a given cell
    static int rowNbr[] = {-1, -1, -1,  0, 0,  1, 1, 1};
    static int colNbr[] = {-1,  0,  1, -1, 1, -1, 0, 1};

    // Mark this cell as visited
    visited.at<uchar>(row, col) = 1;

    // Recur for all connected neighbors
    for (int k = 0; k < 8; ++k)
        if (isSafe(M, row + rowNbr[k], col + colNbr[k], visited, x, matric, count) )
            DFS(M, row + rowNbr[k], col + colNbr[k], visited, x, matric, count);
}

// The main function that returns count of islands in a given boolean
// 2D matrix
int Hls_img::countIslands(Mat M, bool x, Mat matric)
{
    // Make a bool array to mark visited cells.
    // Initially all cells are unvisited
    Mat visited = Mat::zeros(M.rows, M.cols, CV_LOAD_IMAGE_GRAYSCALE);

    // Initialize count as 0 and travel through the all cells of
    // given matrix
    int count = 0;
    for (int i = 0; i < M.rows; ++i)
        for (int j = 0; j < M.cols; ++j)
        {
            Scalar intensity1 = M.at<uchar>(i, j);
            Scalar intensity2 = visited.at<uchar>(i, j);
            if(intensity1.val[0] == x && intensity2.val[0] == 0)
            {
                ++count;
                matric.at<uchar>(i, j) = count;
                DFS(M, i, j, visited, x, matric, count);  // Visit all cells in this island.
            }
        }
    return count;
}

double Hls_img::Distance(int x, int y, int xa, int ya)
{
    double s = 0;
    s = (x - xa)*(x - xa) - (y - ya)*(y - ya);
	return sqrt(s);
}

void Hls_img::Irregularity(Mat image, double irgl0, double irgl1)
{
    irgl0 = irgl1 = 0;
    double max0 = 0, max1 = 0;
    int dem0 = 0, dem1 = 0;
    int sx0 = 0, sy0 = 0, sx1 = 0, sy1 = 0;
    double xa0 =0, ya0 = 0, xa1 = 0, ya1 = 0;
    for (int i=0; i<image.rows; i++)
    {
        for (int j=0; j<image.cols; j++)
        {
            Scalar intensity1 = image.at<uchar>(i, j);
            int x = intensity1.val[0];
            if(x == 0)
            {
                sx0 += i;
                sy0 += j;
                dem0++;
            }
            else
            {
                sx1 += i;
                sy1 += j;
                dem1++;
            }
        }
    }
    xa0 = (double)(sx0/dem0);
    ya0 = (double)(sy0/dem0);
    xa1 = (double)(sx1/dem1);
    ya1 = (double)(sy1/dem1);

    for (int i=0; i<image.rows; i++)
    {
        for (int j=0; j<image.cols; j++)
        {
            Scalar intensity1 = image.at<uchar>(i, j);
            int x = intensity1.val[0];
            if(x == 0)
            {
                if(Distance(i, j, xa0, ya0) > max0)
                    max0 = (double)Hls_img::Distance(i, j, xa0, ya0);
            }
            else
            {
                if(Distance(i, j, xa1, ya1) > max1)
                    max1 = (double)Hls_img::Distance(i, j, xa1, ya1);
            }
        }
    }

    irgl0 = ((1 + sqrt(M_PI)*max0)/(sqrt(dem0))) - 1;
    irgl1 = ((1 + sqrt(M_PI)*max1)/(sqrt(dem1))) - 1;
}

void Hls_img::BienDoiNhiPhan(const char* filename)
{
    //Đọc ảnh màu
    IplImage* img1 = cvLoadImage( filename );

    //Chuyển sang ảnh đa mức xám
    IplImage *img2 = cvCreateImage( cvSize( img1->width, img1->height ), IPL_DEPTH_8U, 1 );

    cvCvtColor( img1, img2, CV_RGB2GRAY );
    IplImage *img3 = cvCreateImage( cvSize( img1->width, img1->height ), IPL_DEPTH_8U, 1 );
    IplImage *img4 = cvCreateImage( cvSize( img1->width, img1->height ), IPL_DEPTH_8U, 1 );

    //Chuyển sang ảnh nhị phân
    cvThreshold(img2, img3, 145, 255, CV_THRESH_BINARY);
    //cvThreshold(img2, img3, 13, 255, CV_THRESH_BINARY);

    //Nguong dong
    cvAdaptiveThreshold(img2, img4, 255, CV_ADAPTIVE_THRESH_MEAN_C, CV_THRESH_BINARY, 5, 5);
    //cvAdaptiveThreshold(img2, img4, 255, CV_ADAPTIVE_THRESH_MEAN_C, CV_THRESH_BINARY, 21, 5);

    Mat im3(img3);
    Mat matric = Mat::zeros(im3.rows, im3.cols, CV_LOAD_IMAGE_GRAYSCALE);
    for(int i = 0; i < im3.rows; i++)
    {
        for(int j = 0; j < im3.cols; j++)
        {
            Scalar intenity = im3.at<uchar>(i, j);
            if(intenity.val[0] == 0)
                matric.at<uchar>(i, j) = 0;
            else
                matric.at<uchar>(i, j) = 1;
        }
    }

    for(int i = 0; i < matric.rows; i++)
    {
        for(int j = 0; j < matric.cols; j++)
        {
            Scalar intenity = matric.at<uchar>(i, j);
            cout<<intenity.val[0]<<"\t";
        }
        cout<<endl;
    }

    Mat matric0 = Mat::zeros(im3.rows, im3.cols, CV_LOAD_IMAGE_GRAYSCALE);
    Mat matric1 = Mat::zeros(im3.rows, im3.cols, CV_LOAD_IMAGE_GRAYSCALE);

    cout<<"Number of islands of 1 is:"<<Hls_img::countIslands(matric, 0, matric0)<<endl;
    for(int i = 0; i < matric0.rows; i++)
    {
        for(int j = 0; j < matric0.cols; j++)
        {
            Scalar intenity = matric0.at<uchar>(i, j);
            cout<<intenity.val[0]<<"\t";
        }
        cout<<endl;
    }

    cout<<"Number of islands of 0 is:"<<countIslands(matric, 1, matric1)<<endl;
    for(int i = 0; i < matric1.rows; i++)
    {
        for(int j = 0; j < matric1.cols; j++)
        {
            Scalar intenity = matric1.at<uchar>(i, j);
            cout<<intenity.val[0]<<"\t";
        }
        cout<<endl;
    }
    namedWindow("1", CV_WINDOW_NORMAL);
    imshow("1", im3);
    waitKey(0);
}


//--------------------------------------------------------------------------
//Entropy, Energy, Inertia for Fourier Transform and Co-occurrence Transform
long double Hls_img::Entropy(Mat image)
{
    long double h = 0;
    long double x = 0;
    for (int i=0; i<image.rows; i++)
    {
        for (int j=0; j<image.cols; j++)
        {
            Scalar intensity1 = image.at<uchar>(i, j);
            x = (long double)intensity1.val[0];
            if(x != 0)
                h += x*log(x);
        }
    }
    return h;
}

long double Hls_img::Energy(Mat image)
{
    long double e = 0;
    long double x = 0;
    for (int i=0; i<image.rows; i++)
    {
        for (int j=0; j<image.cols; j++)
        {
            Scalar intensity1 = image.at<uchar>(i, j);
            x = (long double)intensity1.val[0];
            e += x*x;
        }
    }
    return e;
}

long double Hls_img::Inertia(Mat image)
{
    long double in = 0;
    long double x = 0;
    for (int i=0; i<image.rows; i++)
    {
        for (int j=0; j<image.cols; j++)
        {
            Scalar intensity1 = image.at<uchar>(i, j);
            x = (long double)intensity1.val[0];
            in += (i - j)*(i - j)*x;
        }
    }
    return in;
}


//--------------------------------------------------------------------------
//Co-occurrence Transform
void Hls_img::BienDoiMaTranDongHien(const char * filename)
{
    Mat image = imread(filename, CV_LOAD_IMAGE_GRAYSCALE);
    if( image.empty())
    {
        cout<<"Anh khong ton tai";
        return;
    }
    int row = image.rows, col = image.cols;
    int min_level = 255, max_level = 0;
    for(int i = 0; i < image.rows; i++)
    {
        for(int j = 0; j < image.cols; j++)
        {
            Scalar intensity = image.at<uchar>(i, j);
            int x = intensity.val[0];
            if(x > max_level)
                max_level = x;
            if(x < min_level)
                min_level = x;
        }
    }
    Mat co_matric = Mat::zeros(max_level, max_level, CV_LOAD_IMAGE_GRAYSCALE);
    for(int i = 0; i < image.rows; i++)
    {
        for(int j = 0; j < image.cols; j++)
        {
            if(i - 1 >= 0)
            {
                Scalar intensity1 = image.at<uchar>(i, j);
                Scalar intensity2 = image.at<uchar>(i - 1, j);
                co_matric.at<uchar>(intensity1.val[0], intensity2.val[0])++;
            }
            if(j - 1 >= 0)
            {
                Scalar intensity1 = image.at<uchar>(i, j);
                Scalar intensity2 = image.at<uchar>(i, j - 1);
                co_matric.at<uchar>(intensity1.val[0], intensity2.val[0])++;
            }
            if(j + 1 < col)
            {
                Scalar intensity1 = image.at<uchar>(i, j);
                Scalar intensity2 = image.at<uchar>(i, j + 1);
                co_matric.at<uchar>(intensity1.val[0], intensity2.val[0])++;
            }
            if(i + 1 < row)
            {
                Scalar intensity1 = image.at<uchar>(i, j);
                Scalar intensity2 = image.at<uchar>(i + 1, j + 1);
                co_matric.at<uchar>(intensity1.val[0], intensity2.val[0])++;
            }
        }
    }
    cout<<"He so entropy, h = "<<Hls_img::Entropy(co_matric)<<endl;
    cout<<"He so energy, e = "<<Hls_img::Energy(co_matric)<<endl;
    cout<<"He so inertia, i = "<<Hls_img::Inertia(co_matric)<<endl;

    imshow("gray", image);
    imshow("out", co_matric);
    waitKey(0);
}

void Hls_img::BienDoiMaTranDongHien_8Ke(const char * filename)
{
    Mat image = imread(filename, CV_LOAD_IMAGE_GRAYSCALE);
    if( image.empty())
    {
        cout<<"Anh khong ton tai";
        return;
    }
    int row = image.rows, col = image.cols;
    int min_level = 255, max_level = 0;
    for(int i = 0; i < image.rows; i++)
    {
        for(int j = 0; j < image.cols; j++)
        {
            Scalar intensity = image.at<uchar>(i, j);
            int x = intensity.val[0];
            if(x > max_level)
                max_level = x;
            if(x < min_level)
                min_level = x;
        }
    }
    Mat co_matric = Mat::zeros(max_level, max_level, CV_LOAD_IMAGE_GRAYSCALE);
    for(int i = 0; i < image.rows; i++)
    {
        for(int j = 0; j < image.cols; j++)
        {
            if(i - 1 >= 0)
            {
                if(j - 1 >= 0)
                {
                    Scalar intensity1 = image.at<uchar>(i, j);
                    Scalar intensity2 = image.at<uchar>(i - 1, j - 1);
                    co_matric.at<uchar>(intensity1.val[0], intensity2.val[0])++;
                }
                Scalar intensity1 = image.at<uchar>(i, j);
                Scalar intensity2 = image.at<uchar>(i - 1, j);
                co_matric.at<uchar>(intensity1.val[0], intensity2.val[0])++;
                if(j + 1 >= 0)
                {
                    Scalar intensity1 = image.at<uchar>(i, j);
                    Scalar intensity2 = image.at<uchar>(i - 1, j + 1);
                    co_matric.at<uchar>(intensity1.val[0], intensity2.val[0])++;
                }
            }
            if(j - 1 >= 0)
            {
                Scalar intensity1 = image.at<uchar>(i, j);
                Scalar intensity2 = image.at<uchar>(i, j - 1);
                co_matric.at<uchar>(intensity1.val[0], intensity2.val[0])++;
            }
            if(j + 1 < col)
            {
                Scalar intensity1 = image.at<uchar>(i, j);
                Scalar intensity2 = image.at<uchar>(i, j + 1);
                co_matric.at<uchar>(intensity1.val[0], intensity2.val[0])++;
            }
            if(i + 1 < row)
            {
                if(j - 1 >= 0)
                {
                    Scalar intensity1 = image.at<uchar>(i, j);
                    Scalar intensity2 = image.at<uchar>(i + 1, j - 1);
                    co_matric.at<uchar>(intensity1.val[0], intensity2.val[0])++;
                }
                Scalar intensity1 = image.at<uchar>(i, j);
                Scalar intensity2 = image.at<uchar>(i + 1, j + 1);
                co_matric.at<uchar>(intensity1.val[0], intensity2.val[0])++;
                if(j + 1 >= 0)
                {
                    Scalar intensity1 = image.at<uchar>(i, j);
                    Scalar intensity2 = image.at<uchar>(i + 1, j + 1);
                    co_matric.at<uchar>(intensity1.val[0], intensity2.val[0])++;
                }
            }
        }
    }
    cout<<"He so entropy, h = "<<Hls_img::Entropy(co_matric)<<endl;
    cout<<"He so energy, e = "<<Hls_img::Energy(co_matric)<<endl;
    cout<<"He so inertia, i = "<<Hls_img::Inertia(co_matric)<<endl;

    imshow("gray", image);
    imshow("out", co_matric);
    waitKey(0);
}


//--------------------------------------------------------------------------
//Fourier Transform
void Hls_img::BienDoiFourier(const char* filename)
{
    Mat I = imread(filename, CV_LOAD_IMAGE_GRAYSCALE);
    if( I.empty())
    {
        cout<<"Anh khong ton tai";
        return;
    }

    Mat padded;                            //expand input image to optimal size
    int m = getOptimalDFTSize( I.rows );
    int n = getOptimalDFTSize( I.cols ); // on the border add zero values
    cout<<m<<" x "<<n<<endl;
    copyMakeBorder(I, padded, 0, m - I.rows, 0, n - I.cols, BORDER_CONSTANT, Scalar::all(0));

    Mat planes[] = {Mat_<float>(padded), Mat::zeros(padded.size(), CV_32F)};
    Mat complexI;
    merge(planes, 2, complexI);         // Add to the expanded another plane with zeros

    dft(complexI, complexI);            // this way the result may fit in the source matrix

    // compute the magnitude and switch to logarithmic scale
    // => log(1 + sqrt(Re(DFT(I))^2 + Im(DFT(I))^2))
    split(complexI, planes);                   // planes[0] = Re(DFT(I), planes[1] = Im(DFT(I))
    magnitude(planes[0], planes[1], planes[0]);// planes[0] = magnitude
    Mat magI = planes[0];

    magI += Scalar::all(1);                    // switch to logarithmic scale
    log(magI, magI);

    // crop the spectrum, if it has an odd number of rows or columns
    magI = magI(Rect(0, 0, magI.cols & -2, magI.rows & -2));

    // rearrange the quadrants of Fourier image  so that the origin is at the image center
    int cx = magI.cols/2;
    int cy = magI.rows/2;

    Mat q0(magI, Rect(0, 0, cx, cy));   // Top-Left - Create a ROI per quadrant
    Mat q1(magI, Rect(cx, 0, cx, cy));  // Top-Right
    Mat q2(magI, Rect(0, cy, cx, cy));  // Bottom-Left
    Mat q3(magI, Rect(cx, cy, cx, cy)); // Bottom-Right

    Mat tmp;                           // swap quadrants (Top-Left with Bottom-Right)
    q0.copyTo(tmp);
    q3.copyTo(q0);
    tmp.copyTo(q3);

    q1.copyTo(tmp);                    // swap quadrant (Top-Right with Bottom-Left)
    q2.copyTo(q1);
    tmp.copyTo(q2);

    normalize(magI, magI, 0, 1, CV_MINMAX); // Transform the matrix with float values into a
                                            // viewable image form (float between values 0 and 1).

    cout<<"He so entropy, h = "<<Hls_img::Entropy(magI)<<endl;
    cout<<"He so energy, e = "<<Hls_img::Energy(magI)<<endl;
    cout<<"He so inertia, i = "<<Hls_img::Inertia(magI)<<endl;

    //namedWindow( "Input Image", CV_WINDOW_NORMAL );
    imshow("Input Image", I);    // Show the result
    //namedWindow( "spectrum magnitude", CV_WINDOW_NORMAL );
    imshow("spectrum magnitude", magI);
    waitKey();
}

